Chapter 8 HMSC analysis
8.2 Variance partitioning
# Compute variance partitioning
varpart=computeVariancePartitioning(m)
varpart$vals %>%
as.data.frame() %>%
rownames_to_column(var="variable") %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
group_by(variable) %>%
summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
tt()| variable | mean | sd |
|---|---|---|
| Random: location | 37.900015 | 25.317903 |
| logseqdepth | 56.110626 | 25.796874 |
| sex | 4.937460 | 5.612719 |
| origin | 1.051899 | 1.282563 |
# Basal tree
varpart_tree <- genome_tree
#Varpart table
varpart_table <- varpart$vals %>%
as.data.frame() %>%
rownames_to_column(var="variable") %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(genome=factor(genome, levels=rev(varpart_tree$tip.label))) %>%
mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location"))))
#Phyla
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__"))%>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% varpart_tree$tip.label) %>%
arrange(match(genome, varpart_tree$tip.label)) %>%
mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
column_to_rownames(var = "genome") %>%
dplyr::select(phylum)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__"))%>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% varpart_tree$tip.label) %>%
arrange(match(genome, varpart_tree$tip.label)) %>%
dplyr::select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
dplyr::select(colors) %>%
pull()
# Basal ggtree
varpart_tree <- varpart_tree %>%
force.ultrametric(.,method="extend") %>%
ggtree(., size = 0.3)***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
# Add phylum colors next to the tree tips
varpart_tree <- gheatmap(varpart_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
scale_fill_manual(values=colors_alphabetic)+
labs(fill="Phylum")
#Reset fill scale to use a different colour profile in the heatmap
varpart_tree <- varpart_tree + new_scale_fill()
# Add variance stacked barplot
vertical_tree <- varpart_tree +
scale_fill_manual(values=c("#506a96","#cccccc","#be3e2b","#f6de6c"))+
geom_fruit(
data=varpart_table,
geom=geom_bar,
mapping = aes(x=value, y=genome, fill=variable, group=variable),
pwidth = 2,
offset = 0.05,
width= 1,
orientation="y",
stat="identity")+
labs(fill="Variable")
vertical_tree8.3 Posterior estimates
# Select desired support threshold
support=0.9
negsupport=1-support
# Basal tree
postestimates_tree <- genome_tree
# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
mutate(value = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
pivot_wider(names_from = variable, values_from = value) %>%
#select(genome,sp_vulgaris,area_semi,area_urban,sp_vulgarisxarea_semi,sp_vulgarisxarea_urban,season_spring,season_winter,sp_vulgarisxseason_spring,sp_vulgarisxseason_winter) %>%
column_to_rownames(var="genome")
#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__")) %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% postestimates_tree$tip.label) %>%
arrange(match(genome, postestimates_tree$tip.label)) %>%
mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
column_to_rownames(var = "genome") %>%
dplyr::select(phylum)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__")) %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% postestimates_tree$tip.label) %>%
arrange(match(genome, postestimates_tree$tip.label)) %>%
dplyr::select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
dplyr::select(colors) %>%
pull()
# Basal ggtree
postestimates_tree <- postestimates_tree %>%
force.ultrametric(.,method="extend") %>%
ggtree(., size = 0.3)***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
scale_fill_manual(values=colors_alphabetic)+
labs(fill="Phylum")
#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()
# Add posterior significant heatmap
postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.5, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
labs(fill="Trend")
postestimates_tree +
vexpand(.25, 1) # expand top 8.3.1 Origin
8.3.1.1 Positively associated genomes with domestic cats
getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(trend = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
filter(variable=="originTame", trend=="Positive") %>%
arrange(-value) %>%
left_join(genome_metadata,by=join_by(genome==genome)) %>%
dplyr::select(genome,phylum,class,order,family,species,value) %>%
arrange(phylum, class, family, species)%>%
paged_table()8.3.1.2 Positively associated genomes feral cats
getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(trend = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
filter(variable=="originTame", trend=="Negative") %>%
arrange(value) %>%
left_join(genome_metadata,by=join_by(genome==genome)) %>%
dplyr::select(genome,phylum,class,order,family,species,value) %>%
arrange(phylum,class,family,species)%>%
paged_table()8.3.1.3 Estimate - support plot
estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
filter(variable=="originTame") %>%
pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
dplyr::select(genome,mean)
support <- getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
filter(variable=="originTame") %>%
pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
dplyr::select(genome,support)
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__")) %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% estimate$genome) %>%
arrange(match(genome, estimate$genome)) %>%
dplyr::select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
dplyr::select(colors) %>%
pull()
inner_join(estimate,support,by=join_by(genome==genome)) %>%
mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
mutate(support=ifelse(mean<0,1-support,support)) %>%
left_join(genome_metadata, by = join_by(genome == genome)) %>%
mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
ggplot(aes(x=mean,y=support,color=phylum))+
geom_point(alpha=0.7, shape=16, size=3)+
scale_color_manual(values = phylum_colors) +
geom_vline(xintercept = 0) +
xlim(c(-0.4,0.4)) +
labs(x = "Beta regression coefficient", y = "Posterior probability") +
theme_minimal()+
theme(legend.position = "none")8.3.1.4 Correlations
#Compute the residual correlation matrix
OmegaCor = computeAssociations(m)
# Refernece tree (for sorting genomes)
genome_tree_subset <- genome_tree %>%
keep.tip(., tip=m$spNames)
#Co-occurrence matrix at the animal level
supportLevel = 0.95
toPlot = ((OmegaCor[[1]]$support>supportLevel)
+ (OmegaCor[[1]]$support<(1-supportLevel))>0)*OmegaCor[[1]]$mean
matrix <- toPlot %>%
as.data.frame() %>%
rownames_to_column(var="genome1") %>%
pivot_longer(!genome1, names_to = "genome2", values_to = "cor") %>%
mutate(genome1= factor(genome1, levels=genome_tree_subset$tip.label)) %>%
mutate(genome2= factor(genome2, levels=genome_tree_subset$tip.label)) %>%
ggplot(aes(x = genome1, y = genome2, fill = cor)) +
geom_tile() +
scale_fill_gradient2(low = "#be3e2b",
mid = "#f4f4f4",
high = "#b2b530")+
theme_void()
htree <- genome_tree_subset %>%
force.ultrametric(.,method="extend") %>%
ggtree(.)***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
#create composite figure
grid.arrange(grobs = list(matrix,vtree),
layout_matrix = rbind(c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1)))8.3.1.5 Predict responses (origin)
# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")
gradient = c("domestic","feral")
gradientlength = length(gradient)
#Treatment-specific gradient predictions
pred <- constructGradient(m,
focalVariable = "origin",
non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
predict(m, Gradient = ., expected = TRUE) %>%
do.call(rbind,.) %>%
as.data.frame() %>%
mutate(origin=rep(gradient,1000)) %>%
pivot_longer(!origin,names_to = "genome", values_to = "value")# weights: 9 (4 variable)
initial value 101.072331
final value 91.392443
converged
8.3.1.6 Element level
elements_table <- genome_gifts %>%
to.elements(., GIFT_db) %>%
as.data.frame()
community_elements <- pred %>%
group_by(origin, genome) %>%
mutate(row_id = row_number()) %>%
pivot_wider(names_from = genome, values_from = value) %>%
ungroup() %>%
group_split(row_id) %>%
as.list() %>%
lapply(., FUN = function(x){x %>%
dplyr::select(-row_id) %>%
column_to_rownames(var = "origin") %>%
as.data.frame() %>%
exp() %>%
t() %>%
tss() %>%
to.community(elements_table,.,GIFT_db) %>%
as.data.frame() %>%
rownames_to_column(var="origin")
})
calculate_slope <- function(x) {
lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
coef(lm_fit)[2]
}
element_predictions <- map_dfc(community_elements, function(mat) {
mat %>%
column_to_rownames(var = "origin") %>%
t() %>%
as.data.frame() %>%
rowwise() %>%
mutate(slope = calculate_slope(c_across(everything()))) %>%
dplyr::select(slope) }) %>%
t() %>%
as.data.frame() %>%
set_names(colnames(community_elements[[1]])[-1]) %>%
rownames_to_column(var="iteration") %>%
pivot_longer(!iteration, names_to="trait",values_to="value") %>%
group_by(trait) %>%
summarise(mean=mean(value),
p1 = quantile(value, probs = 0.1),
p9 = quantile(value, probs = 0.9),
positive_support = sum(value > 0)/1000,
negative_support = sum(value < 0)/1000) %>%
arrange(-positive_support)8.3.1.6.1 Positive associated functions at element level
# Positively associated
unique_funct_db<- GIFT_db[c(2,4,5,6)] %>%
distinct(Code_element, .keep_all = TRUE)
element_predictions %>%
filter(mean >0) %>%
arrange(-positive_support) %>%
filter(positive_support>=0.9) %>%
left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
arrange(Domain,Function)%>%
paged_table()#Get community-weighed average GIFTs per sample
genome_counts_row <- genome_counts %>%
mutate_at(vars(-genome),~./sum(.)) %>%
column_to_rownames(., "genome")
#genome_counts_row <- rownames_to_column(genome_counts_row, "genome")
GIFTs_elements_community <- to.community(elements_table,genome_counts_row,GIFT_db)
element_gift_t <- GIFTs_elements_community %>%
t() %>%
# row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
significant_elements<- element_predictions %>%
filter(mean >0) %>%
arrange(-positive_support) %>%
filter(positive_support>=0.9)
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(4,9)], by = join_by(sample == sample))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift_filt %>%
select(-sample)%>%
group_by(Origin) %>%
summarise(across(everything(), mean))%>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element)) Elements Feral Tame Function
1 B0103 0.7566040000 0.7195779000 Nucleic acid biosynthesis_UDP/UTP
2 D0205 0.1081381000 0.1066891000 Polysaccharide degradation_Pectin
3 D0208 0.2168561000 0.2087181000 Polysaccharide degradation_Mixed-Linkage glucans
4 D0504 0.0096225710 0.0100233660 Amino acid degradation_Methionine
5 D0507 0.1643304000 0.1693777000 Amino acid degradation_Leucine
6 D0906 0.0001762601 0.0002054726 Antibiotic degradation_Fosfomycin
8.3.1.6.2 Negatively associated functions at element level
element_predictions %>%
filter(mean <0) %>%
arrange(-negative_support) %>%
filter(negative_support>=0.9) %>%
left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
arrange(Domain,Function)%>%
paged_table()element_gift_t <- GIFTs_elements_community %>%
t() %>%
# row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
significant_elements<- element_predictions %>%
filter(mean <0) %>%
arrange(-negative_support) %>%
filter(negative_support>=0.9)
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(4,9)], by = join_by(sample == sample))
uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift_filt %>%
select(-sample)%>%
group_by(Origin) %>%
summarise(across(everything(), mean))%>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(Elements == Code_element)) Elements Feral Tame Function
1 B0204 0.546283700 0.595230700 Amino acid biosynthesis_Serine
2 B0214 0.211179500 0.236935900 Amino acid biosynthesis_Glutamate
3 B0219 0.011526030 0.010872120 Amino acid biosynthesis_GABA
4 B0302 0.010085548 0.009823828 Amino acid derivative biosynthesis_Betaine
5 B0303 0.233543100 0.236799000 Amino acid derivative biosynthesis_Ectoine
6 B0309 0.076634680 0.071995880 Amino acid derivative biosynthesis_Putrescine
7 B0310 0.055338700 0.046228620 Amino acid derivative biosynthesis_Tryptamine
8 B0401 0.755886800 0.765103900 SCFA biosynthesis_Acetate
9 B0601 0.333605100 0.356787900 Organic anion biosynthesis_Succinate
10 B0603 0.369766100 0.412079900 Organic anion biosynthesis_Citrate
11 B0708 0.440560800 0.464048500 Vitamin biosynthesis_Cobalamin (B12)
12 B0709 0.004493435 0.004851109 Vitamin biosynthesis_Tocopherol/tocotorienol (E)
13 B0804 0.676473200 0.683577900 Aromatic compound biosynthesis_Dipicolinate
14 D0508 0.038700810 0.048971490 Amino acid degradation_Lysine
15 D0512 0.198698300 0.225197200 Amino acid degradation_Histidine
16 D0517 0.039373280 0.046676760 Amino acid degradation_Ornithine
17 D0601 0.063668170 0.064472310 Nitrogen compound degradation_Nitrate
18 D0603 0.009203820 0.010291730 Nitrogen compound degradation_Urate
19 D0606 0.051570030 0.037713990 Nitrogen compound degradation_Allantoin
20 D0610 0.031016560 0.036298430 Nitrogen compound degradation_Methylamine
21 D0611 0.009801617 0.009503118 Nitrogen compound degradation_Phenylethylamine
22 D0612 0.017850890 0.018892270 Nitrogen compound degradation_Hypotaurine
23 D0801 0.002099865 0.001182548 Xenobiotic degradation_Toluene
24 D0802 0.002099865 0.001182548 Xenobiotic degradation_Xylene
25 D0807 0.061055970 0.071058870 Xenobiotic degradation_Catechol
26 D0816 0.091350430 0.092723640 Xenobiotic degradation_Phenylacetate
27 D0817 0.043241270 0.050879070 Xenobiotic degradation_Trans-cinnamate
28 D0903 0.009801617 0.009503118 Antibiotic degradation_Cephalosporin
29 D0908 0.117656700 0.137340600 Antibiotic degradation_Macrolide
positive <- element_predictions %>%
filter(mean >0) %>%
arrange(mean) %>%
filter(positive_support>=0.9) %>%
dplyr::select(-negative_support) %>%
rename(support=positive_support)
negative <- element_predictions %>%
filter(mean <0) %>%
arrange(mean) %>%
filter(negative_support>=0.9) %>%
dplyr::select(-positive_support) %>%
rename(support=negative_support)
bind_rows(positive,negative) %>%
left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
geom_point() +
geom_errorbar() +
xlim(c(-0.04,0.04)) +
geom_vline(xintercept=0) +
scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
theme_minimal() +
labs(x="Regression coefficient",y="Functional trait")table <- bind_rows(positive,negative) %>%
left_join(unique_funct_db,by=join_by(trait==Code_element)) %>%
mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait))))
table %>%
mutate(Element=factor(Element,levels=c(table$Element))) %>%
ggplot(aes(x=mean, y=fct_rev(Element), xmin=p1, xmax=p9, color=Function)) +
geom_point() +
geom_errorbar() +
xlim(c(-0.04,0.04)) +
geom_vline(xintercept=0) +
scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
theme_minimal() +
labs(x="Regression coefficient",y="Functional trait")8.3.1.7 Function level
functions_table <- elements_table %>%
to.functions(., GIFT_db) %>%
as.data.frame()
community_functions <- pred %>%
group_by(origin, genome) %>%
mutate(row_id = row_number()) %>%
pivot_wider(names_from = genome, values_from = value) %>%
ungroup() %>%
group_split(row_id) %>%
as.list() %>%
lapply(., FUN = function(x){x %>%
dplyr::select(-row_id) %>%
column_to_rownames(var = "origin") %>%
as.data.frame() %>%
exp() %>%
t() %>%
tss() %>%
to.community(functions_table,.,GIFT_db) %>%
as.data.frame() %>%
rownames_to_column(var="origin")
})#max-min option
calculate_slope <- function(x) {
lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
coef(lm_fit)[2]
}
function_predictions <- map_dfc(community_functions, function(mat) {
mat %>%
column_to_rownames(var = "origin") %>%
t() %>%
as.data.frame() %>%
rowwise() %>%
mutate(slope = calculate_slope(c_across(everything()))) %>%
dplyr::select(slope) }) %>%
t() %>%
as.data.frame() %>%
set_names(colnames(community_functions[[1]])[-1]) %>%
rownames_to_column(var="iteration") %>%
pivot_longer(!iteration, names_to="trait",values_to="value") %>%
group_by(trait) %>%
summarise(mean=mean(value),
p1 = quantile(value, probs = 0.1),
p9 = quantile(value, probs = 0.9),
positive_support = sum(value > 0)/1000,
negative_support = sum(value < 0)/1000) %>%
arrange(-positive_support)
# Positively associated
function_predictions %>%
filter(mean >0) %>%
arrange(-positive_support) %>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support |
|---|---|---|---|---|---|
| D02 | 0.0081357868 | -0.003637596 | 0.020315293 | 0.825 | 0.175 |
| B08 | 0.0074769039 | -0.003441048 | 0.018287846 | 0.797 | 0.203 |
| B01 | 0.0011491533 | -0.005813419 | 0.007937841 | 0.616 | 0.384 |
| S01 | 0.0008741013 | -0.010932437 | 0.012737002 | 0.563 | 0.437 |
# Negatively associated
function_predictions %>%
filter(mean <0) %>%
arrange(-negative_support) %>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support |
|---|---|---|---|---|---|
| D08 | -1.169835e-03 | -0.0022909469 | -2.422247e-04 | 0.035 | 0.965 |
| B03 | -1.067982e-02 | -0.0178277310 | -3.557808e-03 | 0.045 | 0.955 |
| D06 | -3.274391e-03 | -0.0069404022 | -2.258647e-05 | 0.098 | 0.902 |
| B04 | -8.312913e-03 | -0.0195719001 | 7.637137e-04 | 0.128 | 0.872 |
| B06 | -7.014874e-03 | -0.0169222091 | 2.281875e-03 | 0.169 | 0.831 |
| D07 | -1.194568e-02 | -0.0300927877 | 3.894367e-03 | 0.175 | 0.825 |
| D03 | -4.498865e-03 | -0.0130295373 | 2.733128e-03 | 0.211 | 0.789 |
| D05 | -1.795499e-03 | -0.0075067047 | 3.157093e-03 | 0.223 | 0.777 |
| S03 | -9.537274e-03 | -0.0320097826 | 1.427255e-02 | 0.247 | 0.753 |
| B02 | -3.855277e-03 | -0.0125477755 | 4.170471e-03 | 0.251 | 0.749 |
| S02 | -4.673357e-03 | -0.0146030922 | 2.873079e-03 | 0.304 | 0.696 |
| B07 | -3.854268e-03 | -0.0156450544 | 8.129977e-03 | 0.310 | 0.690 |
| D09 | -1.950214e-03 | -0.0080302387 | 4.401978e-03 | 0.320 | 0.680 |
| B09 | -3.898566e-06 | -0.0005639176 | 5.023253e-04 | 0.366 | 0.634 |
| D01 | -3.475224e-04 | -0.0051398870 | 4.039820e-03 | 0.402 | 0.598 |
| B10 | -2.188595e-05 | -0.0003343497 | 2.311124e-04 | 0.470 | 0.530 |
positive <- function_predictions %>%
filter(mean >0) %>%
arrange(mean) %>%
filter(positive_support>=0.9) %>%
dplyr::select(-negative_support) %>%
rename(support=positive_support)
negative <- function_predictions %>%
filter(mean <0) %>%
arrange(mean) %>%
filter(negative_support>=0.9) %>%
dplyr::select(-positive_support) %>%
rename(support=negative_support)
bind_rows(positive,negative) %>%
left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
geom_point() +
geom_errorbar() +
xlim(c(-0.02,0.02)) +
geom_vline(xintercept=0) +
scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
theme_minimal() +
labs(x="Regression coefficient",y="Functional trait")8.3.2 Sex
8.3.2.1 Negatively associated genomes with male cats
# Compute variance partitioning
varpart=computeVariancePartitioning(m)
varpart$vals %>%
as.data.frame() %>%
rownames_to_column(var="variable") %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
group_by(variable) %>%
summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
tt()| variable | mean | sd |
|---|---|---|
| Random: location | 37.900015 | 25.317903 |
| logseqdepth | 56.110626 | 25.796874 |
| sex | 4.937460 | 5.612719 |
| origin | 1.051899 | 1.282563 |
# Select desired support threshold
support=0.9
negsupport=1-support
# Basal tree
postestimates_tree <- genome_tree
# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
mutate(value = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
pivot_wider(names_from = variable, values_from = value) %>%
column_to_rownames(var="genome")
getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(trend = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
filter(variable=="sexMale", trend=="Negative") %>%
arrange(-value) %>%
left_join(genome_metadata,by=join_by(genome==genome)) %>%
dplyr::select(genome,phylum,class,order,family,species,value) %>%
arrange(phylum, class, family)%>%
tt()| genome | phylum | class | order | family | species | value |
|---|---|---|---|---|---|---|
| bin_Denmark.4 | Actinomycetota | Coriobacteriia | Coriobacteriales | Coriobacteriaceae | Collinsella stercoris | 0.090 |
| bin_Aruba.14 | Actinomycetota | Coriobacteriia | Coriobacteriales | Coriobacteriaceae | Collinsella stercoris | 0.089 |
| bin_Brazil.91 | Bacillota | Bacilli | RF39 | UBA660 | CAG-988 sp003149915 | 0.082 |
| bin_Brazil.76 | Bacillota | Bacilli | RF39 | UBA660 | CAG-877 sp900554305 | 0.017 |
| bin_Malaysia.17 | Bacillota | Bacilli | RF39 | UBA660 | NA | 0.016 |
| bin_Brazil.45 | Bacillota | Bacilli | RF39 | UBA660 | CAG-877 sp900554315 | 0.010 |
| bin_Malaysia.16 | Bacillota_A | Clostridia | Oscillospirales | Acutalibacteraceae | NA | 0.085 |
| bin_Malaysia.78 | Bacillota_A | Clostridia | Oscillospirales | Acutalibacteraceae | Ruminococcus_E bromii_B | 0.081 |
| bin_Brazil.69 | Bacillota_A | Clostridia | Oscillospirales | Acutalibacteraceae | Clostridium_A leptum | 0.020 |
| bin_Brazil.83 | Bacillota_A | Clostridia | Lachnospirales | Anaerotignaceae | NA | 0.069 |
| bin_Denmark.63 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Anaerostipes hadrus | 0.091 |
| bin_Spain.11 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.084 |
| bin_Malaysia.3 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.077 |
| bin_Malaysia.45 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Blautia_A wexlerae | 0.077 |
| bin_Brazil.1 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.071 |
| bin_Malaysia.7 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | CAG-317 sp900543415 | 0.070 |
| bin_Denmark.118 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Dorea_A longicatena | 0.066 |
| bin_Brazil.3 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Faecalimonas sp900555395 | 0.065 |
| bin_Malaysia.52 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Fusicatenibacter saccharivorans | 0.065 |
| bin_Brazil.89 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Faecalimonas sp900550235 | 0.062 |
| bin_Denmark.34 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Fusicatenibacter saccharivorans | 0.062 |
| bin_Malaysia.73 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.062 |
| bin_Brazil.166 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.057 |
| bin_Brazil.165 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Blautia_A caecimuris | 0.055 |
| bin_Denmark.19 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Faecalimonas sp900556835 | 0.053 |
| bin_Spain.53 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.053 |
| bin_Brazil.54 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Dorea_B phocaeensis | 0.051 |
| bin_Brazil.157 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.046 |
| bin_Brazil.32 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Mediterraneibacter torques | 0.042 |
| bin_Brazil.8 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Blautia_A sp900547615 | 0.042 |
| bin_Spain.6 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Schaedlerella glycyrrhizinilytica | 0.039 |
| bin_Denmark.109 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Eisenbergiella sp900550285 | 0.034 |
| bin_Brazil.113 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Anaerobutyricum sp900754855 | 0.031 |
| bin_Brazil.17 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.030 |
| bin_Spain.37 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | UMGS1472 sp900552095 | 0.029 |
| bin_Denmark.52 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Ruminococcus_B gnavus | 0.023 |
| bin_Brazil.56 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Ruminococcus_A sp000432335 | 0.022 |
| bin_Spain.67 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Ruminococcus_B sp900544395 | 0.021 |
| bin_Brazil.28 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.019 |
| bin_Brazil.93 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.019 |
| bin_Malaysia.110 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Catenibacillus sp900557175 | 0.018 |
| bin_Aruba.43 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Roseburia sp900548205 | 0.017 |
| bin_Malaysia.108 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.017 |
| bin_Brazil.116 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.016 |
| bin_Malaysia.86 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.015 |
| bin_Brazil.99 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Lachnospira sp900552795 | 0.013 |
| bin_CaboVerde.61 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Roseburia sp900548205 | 0.013 |
| bin_Brazil.97 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | UBA9502 sp900538475 | 0.010 |
| bin_Malaysia.46 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Blautia sp003287895 | 0.010 |
| bin_Brazil.50 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Lachnospira sp003451515 | 0.009 |
| bin_Denmark.66 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Blautia sp900120295 | 0.009 |
| bin_Spain.60 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.009 |
| bin_Malaysia.125 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Blautia sp900539145 | 0.008 |
| bin_CaboVerde.34 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.005 |
| bin_Aruba.13 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | CAG-81 sp000435795 | 0.002 |
| bin_Brazil.63 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.002 |
| bin_Denmark.62 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Clostridium_Q sp000435655 | 0.002 |
| bin_Brazil.105 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.001 |
| bin_Denmark.20 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Enterocloster sp001517625 | 0.001 |
| bin_Malaysia.98 | Bacillota_A | Clostridia | Oscillospirales | Oscillospiraceae | NA | 0.066 |
| bin_Malaysia.9 | Bacillota_A | Clostridia | Oscillospirales | Oscillospiraceae | Dysosmobacter welbionis | 0.065 |
| bin_Aruba.28 | Bacillota_A | Clostridia | Oscillospirales | Oscillospiraceae | NA | 0.044 |
| bin_Brazil.177 | Bacillota_A | Clostridia | Oscillospirales | Oscillospiraceae | NA | 0.044 |
| bin_Malaysia.116 | Bacillota_A | Clostridia | Oscillospirales | Oscillospiraceae | Flavonifractor plautii | 0.043 |
| bin_Malaysia.34 | Bacillota_A | Clostridia | Oscillospirales | Oscillospiraceae | Lawsonibacter sp000177015 | 0.040 |
| bin_Brazil.137 | Bacillota_A | Clostridia | Oscillospirales | Oscillospiraceae | NA | 0.027 |
| bin_Aruba.52 | Bacillota_A | Clostridia | Tissierellales | Peptoniphilaceae | NA | 0.094 |
| bin_Aruba.31 | Bacillota_A | Clostridia | Oscillospirales | Ruminococcaceae | Negativibacillus sp000435195 | 0.080 |
| bin_Malaysia.102 | Bacillota_A | Clostridia | Oscillospirales | Ruminococcaceae | NA | 0.080 |
| bin_Malaysia.30 | Bacillota_A | Clostridia | Oscillospirales | Ruminococcaceae | Phocea massiliensis | 0.055 |
| bin_Denmark.72 | Bacillota_C | Negativicutes | Selenomonadales | Selenomonadaceae | Megamonas funiformis | 0.082 |
| bin_Denmark.85 | Bacillota_C | Negativicutes | Selenomonadales | Selenomonadaceae | NA | 0.077 |
| bin_Brazil.5 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotella sp900540415 | 0.080 |
| bin_Malaysia.18 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotellamassilia sp000437675 | 0.051 |
| bin_Brazil.48 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Bacteroides sp900766005 | 0.050 |
| bin_Malaysia.117 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotellamassilia sp900541335 | 0.047 |
| bin_CaboVerde.18 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotella copri | 0.038 |
| bin_Aruba.10 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotella sp900544825 | 0.037 |
| bin_Malaysia.64 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Paraprevotella clara | 0.033 |
| bin_Brazil.163 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotella lascolaii | 0.029 |
| bin_Spain.4 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Phocaeicola sp000436795 | 0.027 |
| bin_Spain.48 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Bacteroides fragilis | 0.019 |
| bin_Brazil.96 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Bacteroides uniformis | 0.014 |
| bin_Malaysia.105 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Phocaeicola massiliensis | 0.014 |
| bin_Spain.21 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Phocaeicola coprophilus | 0.012 |
| bin_Malaysia.121 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Bacteroides stercoris | 0.011 |
| bin_Brazil.43 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Phocaeicola vulgatus | 0.010 |
| bin_Denmark.30 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Phocaeicola sp900546645 | 0.010 |
| bin_Brazil.103 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Phocaeicola plebeius_A | 0.009 |
| bin_Malaysia.77 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Phocaeicola sp900542985 | 0.006 |
| bin_Brazil.6 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Phocaeicola coprocola | 0.005 |
| bin_Brazil.110 | Bacteroidota | Bacteroidia | Bacteroidales | Barnesiellaceae | Barnesiella intestinihominis | 0.024 |
| bin_Malaysia.13 | Bacteroidota | Bacteroidia | Bacteroidales | Muribaculaceae | NA | 0.075 |
| bin_Spain.23 | Bacteroidota | Bacteroidia | Bacteroidales | Muribaculaceae | CAG-279 sp900541935 | 0.051 |
| bin_CaboVerde.37 | Bacteroidota | Bacteroidia | Bacteroidales | Porphyromonadaceae | NA | 0.058 |
| bin_Brazil.11 | Bacteroidota | Bacteroidia | Bacteroidales | Rikenellaceae | Alistipes putredinis | 0.045 |
| bin_Malaysia.131 | Bacteroidota | Bacteroidia | Bacteroidales | Rikenellaceae | Alistipes putredinis | 0.043 |
| bin_Brazil.38 | Bacteroidota | Bacteroidia | Bacteroidales | Rikenellaceae | Alistipes communis | 0.039 |
| bin_Brazil.86 | Bacteroidota | Bacteroidia | Bacteroidales | Rikenellaceae | Alistipes dispar | 0.030 |
| bin_Brazil.138 | Bacteroidota | Bacteroidia | Bacteroidales | Tannerellaceae | Parabacteroides sp000436495 | 0.023 |
| bin_Brazil.124 | Bacteroidota | Bacteroidia | Bacteroidales | Tannerellaceae | NA | 0.019 |
| bin_Brazil.160 | Bacteroidota | Bacteroidia | Bacteroidales | Tannerellaceae | Parabacteroides distasonis | 0.013 |
| bin_CaboVerde.10 | Campylobacterota | Campylobacteria | Campylobacterales | Campylobacteraceae | NA | 0.070 |
| bin_Brazil.68 | Campylobacterota | Campylobacteria | Campylobacterales | Campylobacteraceae | Campylobacter_D helveticus | 0.049 |
| bin_Spain.26 | Campylobacterota | Campylobacteria | Campylobacterales | Campylobacteraceae | Campylobacter_D upsaliensis | 0.047 |
| bin_Brazil.145 | Cyanobacteria | Vampirovibrionia | Gastranaerophilales | Gastranaerophilaceae | NA | 0.044 |
| bin_Brazil.140 | Fusobacteriota | Fusobacteriia | Fusobacteriales | Fusobacteriaceae | Fusobacterium_A sp900543175 | 0.029 |
| bin_Malaysia.6 | Fusobacteriota | Fusobacteriia | Fusobacteriales | Fusobacteriaceae | Fusobacterium_B sp900545035 | 0.015 |
| bin_Brazil.159 | Fusobacteriota | Fusobacteriia | Fusobacteriales | Fusobacteriaceae | Fusobacterium_B sp900541465 | 0.012 |
| bin_Brazil.146 | Pseudomonadota | Alphaproteobacteria | RF32 | CAG-239 | CAG-495 sp001917125 | 0.074 |
| bin_Malaysia.50 | Pseudomonadota | Gammaproteobacteria | Enterobacterales | Succinivibrionaceae | NA | 0.084 |
| bin_Brazil.82 | Pseudomonadota | Gammaproteobacteria | Enterobacterales | Succinivibrionaceae | Anaerobiospirillum succiniciproducens | 0.046 |
| bin_CaboVerde.33 | Pseudomonadota | Gammaproteobacteria | Enterobacterales | Succinivibrionaceae | Anaerobiospirillum sp900543125 | 0.045 |
| bin_CaboVerde.55 | Pseudomonadota | Gammaproteobacteria | Enterobacterales | Succinivibrionaceae | Anaerobiospirillum_A thomasii | 0.036 |
| bin_Brazil.111 | Pseudomonadota | Gammaproteobacteria | Enterobacterales | Succinivibrionaceae | Succinivibrio sp000431835 | 0.024 |
8.3.2.2 Estimate - support plot
estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
filter(variable=="sexMale") %>%
pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
dplyr::select(genome,mean)
support <- getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
filter(variable=="sexMale") %>%
pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
dplyr::select(genome,support)
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__")) %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% estimate$genome) %>%
arrange(match(genome, estimate$genome)) %>%
dplyr::select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
dplyr::select(colors) %>%
pull()
inner_join(estimate,support,by=join_by(genome==genome)) %>%
mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
mutate(support=ifelse(mean<0,1-support,support)) %>%
left_join(genome_metadata, by = join_by(genome == genome)) %>%
mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
ggplot(aes(x=mean,y=support,color=phylum))+
geom_point(alpha=0.7, shape=16, size=3)+
scale_color_manual(values = phylum_colors) +
geom_vline(xintercept = 0) +
xlim(c(-0.4,0.4)) +
labs(x = "Beta regression coefficient", y = "Posterior probability") +
theme_minimal()+
theme(legend.position = "none")8.3.2.3 Predict responses
# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")
gradient = c("Male","Female","Unknown")
gradientlength = length(gradient)
#Treatment-specific gradient predictions
pred <- constructGradient(m,
focalVariable = "sex",
non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
predict(m, Gradient = ., expected = TRUE) %>%
do.call(rbind,.) %>%
as.data.frame() %>%
mutate(sex=rep(gradient,1000)) %>%
pivot_longer(!sex,names_to = "genome", values_to = "value")# weights: 4 (3 variable)
initial value 63.769541
final value 61.728471
converged
8.3.2.4 Element level
elements_table <- genome_gifts %>%
to.elements(., GIFT_db) %>%
as.data.frame()
community_elements <- pred %>%
group_by(sex, genome) %>%
mutate(row_id = row_number()) %>%
pivot_wider(names_from = genome, values_from = value) %>%
ungroup() %>%
group_split(row_id) %>%
as.list() %>%
lapply(., FUN = function(x){x %>%
dplyr::select(-row_id) %>%
column_to_rownames(var = "sex") %>%
as.data.frame() %>%
exp() %>%
t() %>%
tss() %>%
to.community(elements_table,.,GIFT_db) %>%
as.data.frame() %>%
rownames_to_column(var="sex")
})
calculate_slope <- function(x) {
lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
coef(lm_fit)[2]
}
element_predictions <- map_dfc(community_elements, function(mat) {
mat %>%
column_to_rownames(var = "sex") %>%
t() %>%
as.data.frame() %>%
rowwise() %>%
mutate(slope = calculate_slope(c_across(everything()))) %>%
dplyr::select(slope) }) %>%
t() %>%
as.data.frame() %>%
set_names(colnames(community_elements[[1]])[-1]) %>%
rownames_to_column(var="iteration") %>%
pivot_longer(!iteration, names_to="trait",values_to="value") %>%
group_by(trait) %>%
summarise(mean=mean(value),
p1 = quantile(value, probs = 0.1),
p9 = quantile(value, probs = 0.9),
positive_support = sum(value > 0)/1000,
negative_support = sum(value < 0)/1000) %>%
arrange(-positive_support)8.3.2.4.1 Positive associated functions at element level
# Positively associated
unique_funct_db<- GIFT_db[c(2,4,5,6)] %>%
distinct(Code_element, .keep_all = TRUE)
element_predictions %>%
filter(mean >0) %>%
arrange(-positive_support) %>%
filter(positive_support>=0.9) %>%
left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
arrange(Domain,Function)# A tibble: 27 × 9
trait mean p1 p9 positive_support negative_support Domain Function Element
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr>
1 B0212 0.0295 0.00137 0.0554 0.909 0.091 Biosynthesis Amino acid biosynthesis Arginine
2 B0220 0.00838 0.000520 0.0159 0.909 0.091 Biosynthesis Amino acid biosynthesis Beta-alanine
3 B0310 0.0189 0.00558 0.0325 0.959 0.041 Biosynthesis Amino acid derivative biosynthesis Tryptamine
4 B0303 0.0146 0.000861 0.0264 0.913 0.087 Biosynthesis Amino acid derivative biosynthesis Ectoine
5 B0307 0.0408 0.00190 0.0709 0.906 0.094 Biosynthesis Amino acid derivative biosynthesis Spermidine
6 B0901 0.00110 0.000165 0.00258 0.917 0.083 Biosynthesis Metallophore biosynthesis Staphyloferrin
7 B0104 0.0328 0.00473 0.0561 0.925 0.075 Biosynthesis Nucleic acid biosynthesis CDP/CTP
8 B0105 0.0245 0.00614 0.0419 0.924 0.076 Biosynthesis Nucleic acid biosynthesis ADP/ATP
9 B0106 0.0139 0.000246 0.0265 0.901 0.099 Biosynthesis Nucleic acid biosynthesis GDP/GTP
10 B0605 0.0618 0.0328 0.0910 0.966 0.034 Biosynthesis Organic anion biosynthesis D-lactate
# ℹ 17 more rows
8.3.2.4.2 Negatively associated functions at element level
element_predictions %>%
filter(mean <0) %>%
arrange(-negative_support) %>%
filter(negative_support>=0.9) %>%
left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
arrange(Domain,Function)# A tibble: 19 × 9
trait mean p1 p9 positive_support negative_support Domain Function Element
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr> <chr>
1 B0215 -0.0600 -0.0910 -0.0281 0.028 0.972 Biosynthesis Amino acid biosynthesis Histidine
2 B0216 -0.0477 -0.0710 -0.0229 0.031 0.969 Biosynthesis Amino acid biosynthesis Tryptophan
3 B0208 -0.0286 -0.0528 -0.000896 0.093 0.907 Biosynthesis Amino acid biosynthesis Valine
4 B0209 -0.0274 -0.0539 -0.0000683 0.1 0.9 Biosynthesis Amino acid biosynthesis Isoleucine
5 B1012 -0.00822 -0.0129 -0.00371 0.007 0.993 Biosynthesis Antibiotic biosynthesis Fosfomycin
6 B0704 -0.0585 -0.0919 -0.0262 0.019 0.981 Biosynthesis Vitamin biosynthesis Pantothenate (B5)
7 B0711 -0.0331 -0.0528 -0.0143 0.041 0.959 Biosynthesis Vitamin biosynthesis Menaquinone (K2)
8 B0710 -0.0177 -0.0297 -0.00521 0.059 0.941 Biosynthesis Vitamin biosynthesis Phylloquinone (K1)
9 B0703 -0.0134 -0.0263 -0.00142 0.077 0.923 Biosynthesis Vitamin biosynthesis Niacin (B3)
10 D0516 -0.0000411 -0.0000755 -0.00000151 0.008 0.992 Degradation Amino acid degradation Beta-alanine
11 D0503 -0.000644 -0.00129 -0.0000797 0.056 0.944 Degradation Amino acid degradation Cysteine
12 D0906 -0.00350 -0.00696 -0.000488 0.04 0.96 Degradation Antibiotic degradation Fosfomycin
13 D0613 -0.0604 -0.0970 -0.0233 0.038 0.962 Degradation Nitrogen compound degradation Taurine
14 D0209 -0.0461 -0.0729 -0.0201 0.017 0.983 Degradation Polysaccharide degradation Xylans
15 D0205 -0.0200 -0.0312 -0.00907 0.02 0.98 Degradation Polysaccharide degradation Pectin
16 D0210 -0.0160 -0.0316 -0.00271 0.056 0.944 Degradation Polysaccharide degradation Beta-mannan
17 D0202 -0.0170 -0.0300 -0.00368 0.064 0.936 Degradation Polysaccharide degradation Xyloglucan
18 D0201 -0.0194 -0.0323 -0.00546 0.069 0.931 Degradation Polysaccharide degradation Cellulose
19 D0203 -0.0154 -0.0284 -0.00110 0.093 0.907 Degradation Polysaccharide degradation Starch
positive <- element_predictions %>%
filter(mean >0) %>%
arrange(mean) %>%
filter(positive_support>=0.9) %>%
dplyr::select(-negative_support) %>%
rename(support=positive_support)
negative <- element_predictions %>%
filter(mean <0) %>%
arrange(mean) %>%
filter(negative_support>=0.9) %>%
dplyr::select(-positive_support) %>%
rename(support=negative_support)
bind_rows(positive,negative) %>%
left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
geom_point() +
geom_errorbar() +
# xlim(c(-0.04,0.04)) +
geom_vline(xintercept=0) +
# scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
theme_minimal() +
labs(x="Regression coefficient",y="Functional trait")8.3.2.5 Function level
functions_table <- elements_table %>%
to.functions(., GIFT_db) %>%
as.data.frame()
community_functions <- pred %>%
group_by(sex, genome) %>%
mutate(row_id = row_number()) %>%
pivot_wider(names_from = genome, values_from = value) %>%
ungroup() %>%
group_split(row_id) %>%
as.list() %>%
lapply(., FUN = function(x){x %>%
dplyr::select(-row_id) %>%
column_to_rownames(var = "sex") %>%
as.data.frame() %>%
exp() %>%
t() %>%
tss() %>%
to.community(functions_table,.,GIFT_db) %>%
as.data.frame() %>%
rownames_to_column(var="sex")
})#max-min option
calculate_slope <- function(x) {
lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
coef(lm_fit)[2]
}
function_predictions <- map_dfc(community_functions, function(mat) {
mat %>%
column_to_rownames(var = "sex") %>%
t() %>%
as.data.frame() %>%
rowwise() %>%
mutate(slope = calculate_slope(c_across(everything()))) %>%
dplyr::select(slope) }) %>%
t() %>%
as.data.frame() %>%
set_names(colnames(community_functions[[1]])[-1]) %>%
rownames_to_column(var="iteration") %>%
pivot_longer(!iteration, names_to="trait",values_to="value") %>%
group_by(trait) %>%
summarise(mean=mean(value),
p1 = quantile(value, probs = 0.1),
p9 = quantile(value, probs = 0.9),
positive_support = sum(value > 0)/1000,
negative_support = sum(value < 0)/1000) %>%
arrange(-positive_support)
# Positively associated
function_predictions %>%
filter(mean >0) %>%
arrange(-positive_support) %>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support |
|---|---|---|---|---|---|
| D07 | 0.0370879175 | 1.865579e-02 | 0.054600433 | 0.990 | 0.010 |
| B03 | 0.0132361030 | 2.560580e-03 | 0.022378474 | 0.924 | 0.076 |
| B01 | 0.0112990801 | 1.536807e-03 | 0.021070397 | 0.922 | 0.078 |
| B09 | 0.0008303995 | 7.280411e-06 | 0.001905481 | 0.901 | 0.099 |
| D06 | 0.0042811880 | -2.679345e-05 | 0.009422817 | 0.898 | 0.102 |
| D03 | 0.0081369808 | -1.021315e-03 | 0.017939832 | 0.879 | 0.121 |
| B04 | 0.0089402186 | -1.852843e-03 | 0.019829339 | 0.870 | 0.130 |
| D09 | 0.0045467180 | -4.073236e-03 | 0.012796230 | 0.819 | 0.181 |
| D05 | 0.0033674482 | -7.647938e-03 | 0.011635638 | 0.809 | 0.191 |
| B06 | 0.0097564043 | -5.525707e-03 | 0.025204295 | 0.807 | 0.193 |
| S03 | 0.0152479121 | -1.896561e-02 | 0.043166598 | 0.799 | 0.201 |
| D08 | 0.0003058347 | -3.947465e-04 | 0.001036747 | 0.685 | 0.315 |
| D01 | 0.0004439636 | -5.283954e-03 | 0.006055686 | 0.646 | 0.354 |
# Negatively associated
function_predictions %>%
filter(mean <0) %>%
arrange(-negative_support) %>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support |
|---|---|---|---|---|---|
| B10 | -0.0006188166 | -0.0009435699 | -0.0002871661 | 0.029 | 0.971 |
| D02 | -0.0148823215 | -0.0276192261 | -0.0021080715 | 0.074 | 0.926 |
| B07 | -0.0093042984 | -0.0214646905 | 0.0021940396 | 0.134 | 0.866 |
| S01 | -0.0128671477 | -0.0284627187 | 0.0050748248 | 0.139 | 0.861 |
| S02 | -0.0033662258 | -0.0116276565 | 0.0056703820 | 0.179 | 0.821 |
| B02 | -0.0035628916 | -0.0153096303 | 0.0077135414 | 0.365 | 0.635 |
| B08 | -0.0024022736 | -0.0160211150 | 0.0099821510 | 0.514 | 0.486 |
positive <- function_predictions %>%
filter(mean >0) %>%
arrange(mean) %>%
filter(positive_support>=0.9) %>%
dplyr::select(-negative_support) %>%
rename(support=positive_support)
negative <- function_predictions %>%
filter(mean <0) %>%
arrange(mean) %>%
filter(negative_support>=0.9) %>%
dplyr::select(-positive_support) %>%
rename(support=negative_support)
bind_rows(positive,negative) %>%
left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
geom_point() +
geom_errorbar() +
# xlim(c(-0.02,0.02)) +
geom_vline(xintercept=0) +
scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
theme_minimal() +
labs(x="Regression coefficient",y="Functional trait")